clear all;

cd ..; cd ..; cd ..;
load('code_main_model\6_implied_w_r\computation_results\data.mat','gb_ini','gtheta_ini','V_ini',...
    'gridtheta','gridb','gridz','nmkt','nz','nb','ntheta','kappa_ini');
cd 'code_figure_table\OA\Figure OA8_9';

kappa=kappa_ini;
V=squeeze(V_ini(1,:,:,:,:));
gb=gb_ini;
gtheta=gtheta_ini;

mig_1to2(1:nz,1:nb,1:ntheta)=0;
mig_2to1(1:nz,1:nb,1:ntheta)=0;

i_mkt=1;
i_mkt_to=2;
for i_z=1:nz
    for i_b=1:nb
        for i_theta=1:ntheta
            b=gb(i_mkt, i_z, i_b, i_theta);
            theta=gtheta(i_mkt, i_z, i_b, i_theta);
            m=b*theta;
            b_M=b-kappa;
            m_M=m-kappa;
            if (m_M<0)
                continue;
            else
                if b_M<1e-10
                    theta_M=1;
                else
                    theta_M=m_M/b_M;
                end
                
                V_stay=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(V(i_mkt,i_z,:,:)), b, theta);
                V_mig=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(V(i_mkt_to,i_z,:,:)), b_M, theta_M);

                if V_mig>V_stay
                    mig_1to2(i_z, i_b, i_theta)=1;
                end
            end
           
        end
    end
end


i_mkt=2;
i_mkt_to=1;
for i_z=1:nz
    for i_b=1:nb
        for i_theta=1:ntheta
            b=gb(i_mkt, i_z, i_b, i_theta);
            theta=gtheta(i_mkt, i_z, i_b, i_theta);
            m=b*theta;
            b_M=b-kappa;
            m_M=m-kappa;
            if (m_M<0)
                continue;
            else
                if b_M<1e-10
                    theta_M=1;
                else
                    theta_M=m_M/b_M;
                end
                
                V_stay=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(V(i_mkt,i_z,:,:)), b, theta);
                V_mig=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(V(i_mkt_to,i_z,:,:)), b_M, theta_M);

                if V_mig>V_stay
                    mig_2to1(i_z, i_b, i_theta)=1;
                end
            end
           
        end
    end
end

save('mig.mat','mig_1to2','mig_2to1','gridtheta','gridb','gridz','nmkt','nz','nb','ntheta','kappa');